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In this article, we compare the results ol non-equilibrium (NEMD) and equilibrium (EMD) molec- 
ular dynamics methods to compute the thermal conductance at the interlace between solids. We 
propose to probe the thermal conductance using equilibrium simulations measuring the decay ol 
the thermally induced energy fluctuations ol each solid. We also show that NEMD and EMD 
give generally speaking inconsistent results lor the thermal conductance: Green Kubo simulations 
probe the Landauer conductance between two solids which assumes phonons on both sides ol the 
interlace to be at equilibrium. On the other hand, we show that NEMD give access to the out-of- 
equilibrium interlacial conductance consistent with the interlacial flux describing phonon transport 
in each solid. The difference may be large and reaches typically a lactor 5 lor interlaces between 
usual semi-conductors. We analyze finite size effects lor the two determinations ol the interlacial 
thermal conductance, and show that the equilibrium simulations suffer Irom severe size effects as 
compared to NEMD. We also compare the predictions ol the two above mentioned methods -EMD 
and NEMD- regarding the interlacial conductance ol a series ol mass mismatched Lennard- Jones 
solids. We show that the Kapitza conductance obtained with EMD can be well described using 
the classical diffuse mismatch model (DMM). On the other hand, NEMD simulations results are 
consistent with a out-ol-equilibrium generalisation ol the acoustic mismatch model (AMM). These 
considerations are important in rationalizing previous results obtained using molecular dynamics, 
and help in pinpointing the physical scattering mechanisms taking place at atomically perfect in- 
terlaces between solids, which is a prerequesite to understand interlacial heat transler across real 
interlaces. 
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I. INTRODUCTION 

Kapitza conductance controls heat transfer at submi- 
cronic length scales in heterogeneous and nanostructured 
materials. For instance in superlattices, which are made 
of an arrangement of alternating solid layers, the Kapitza 
conductance at the interface between the solids controls 
the overall conductivity of the superlattice when the in- 
ternal conductance of the solid layers is large — . Un- 
derstanding the value of the Kapitza conductance at the 
interface between solids may thus help in defining direc- 
tions to minimize or on the contrary maximise the con- 
ductivity of the superlattice, with respective applications 
in energy conversion devices and thermal management. 
During the last decade, ultrafast measurements tech- 
niques have been developed so that the Kapitza 
conductance at the interface between a number of 
metal/dielectrics and dielectrics/dielectrics solids has 
been characterized 3 - - — . Similarly, ultrafast LASER spec- 
troscopy may also be used to measure the Kapitza con- 
ductance between a metal and a solid matrix which can 
be amorphous^. All the above-mentioned experiments 
have concluded that the Kapitza conductance is poorly 
described by the classical AMM and DMM models, with 
sometimes a difference reaching an order of magnitude. 
Also the temperature dependence predicted by the clas- 



sical models is wrong, with experiments and simulations 
pointing at a linear increase of the conductance with the 
temperature^£&i2 when the theories predict a constant 
value at least if interfacial scattering is supposed to be 
elastic. These discrepancies may be partly explained by 
the state of the interface between real materials whose 
imperfections may enhance inelastic scattering, thus cre- 
ating additional energy channels compared with the sit- 
uation of an ideal interface. In this context theoretical 
modeling may help in pinpointing the physical relevant 
mechanisms ruling heat transfer across ideal interfaces. 
To this end, different techniques have been employed in- 
cluding lattice dynamics^, Green Function^ and molec- 
ular dynamics (MD) 11 ' 12 . The latter is a promising 
method as it is relatively easy to use and it makes no as- 
sumption regarding interfacial heat tranport except the 
classical nature of the energy carriers, a reasonable as- 
sumption close to the Debye temperature of the softer 
solid. However, even for perfect interfaces no agreement 
has been found between the MD results and the classical 
AMM and DMM models!^. As bulk transport coeffi- 
cients, two routes may be followed to determine the in- 
terfacial conductance between classical solids : either the 
system is driven out-of-equilibrium by creating an inter- 
facial flux using two heat reservoirs on both sides of the 
interfac o 12 ! 13 or the kinetics of thermally induced flue- 
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tuations of the interfacial flux may be recorded around 
the equilibrium situation where the two solids are at the 
same temperature^—. This latter method relies on the 
generalisation of the Green-Kubo formulae to interfacial 
transport coefficients^. Contrary to the case of the ther- 
mal conductivity, no agreement has been found between 
these two methods even when considering simple systems 
such as the interface between Lennard-Jones solids^. 
In this article, we propose a new method to determine the 
interfacial conductance in the spirit of the EMD method 
using the energy autocorrelation function of each solid on 
both sides of the interface. This may solve practical prob- 
lems frequently encountered in equilibrium simulations 
when a plateau in the integral of the relevant correlation 
function should be identified, which often leads to prac- 
tical difficulties. We explain the discrepancies between 
the EMD and the NEMD simulations determination of 
the interfacial heat conductance. We show that the EMD 
yields the Landauer conductance which assumes phonons 
on both sides of the interface to have equilibrium distri- 
bution. On the other hand, we will show that the con- 
ductance measured in NEMD is well described by the 
general expression of Simons which accounts for the out- 
of-equilibrium phonon distribution consistent with the 
created heat flu x 17 ' 27 . Thus we conclude that the two 
methods give intrinsically inconsistent values of the inter- 
facial conductance. The difference is important and may 
reach nearly an order of magnitude for solids displaying 
moderate acoustic mismatch. We analyze also the finite 
size effects in the two methods and show that EMD suf- 
fers from stronger size effects than NEMD. Finally, we 
analyze the interfacial conductance at the interface be- 
tween a series of mass-mismatched Lennard-Jones solids 
using both methods. We show that the classical DMM 
model provides a good description of the EMD data. On 
the other hand, both the AMM and DMM models fail 
to predict the conductance obtained in NEMD. A good 
agreement is found if we extend the AMM model by ac- 
counting for the out-of-equilibrium phonon distribution 
consistent with the imposed interfacial flux. 
The article is structured as follows: in the section [ill we 
first review the basics of interfacial heat transport; We 
discuss the difference between the Landauer conductance 
which assumes the energy carriers to be described locally 
by equilibrium distribution functions and the general ex- 
pression proposed by Simons. In the section Hill we show 
the connection between the Landauer conductance and 
the decay of the energy autocorrelation function in each 
solid. This allows us to propose an alternate expression 
to measure the interfacial conductance using EMD. This 
methodology is applied in the section IIVI where we an- 
alyze the case of the interface between Lennard-Jones 
solids having a variable mass contrast. We also compare 
the conductance obtained using the two methods with 
the different theoretical predictions discussed in the sec- 
tion |nl We discuss the consequences of this work in the 
Conclusion. 
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FIG. 1. (Color online) Temperature jump across the interface 
between two media crossed by an interfacial flux q. The tem- 
perature profile in each medium is schematically represented 
by red solid lines. 

II. THEORY 

In this section, we briefly review the basic definitions 
of the interfacial conductance and we discuss its relation 
with the phonon distribution on both sides of the 
interface. The equations derived in this section are not 
completely new but they are reviewed for the sake of 
completeness. In particular, we review the expression 
first proposed by Simons^ which accounts for the 
out-of-equilibrium phonon distribution consistent with 
the interfacial flux. 

Consider the interface between two media 1 and 2 as 
sketched in fig. [TJ The interfacial conductance G between 
these two media is defined in terms of the ratio 

G = q/(T 2 -T 1 ) (1) 

where q is the heat flux flowing across the interface from 
medium 2 to 1, and Ti denotes the temperature of the 
medium i in the vicinity of the interface. The interfacial 
conductance eq. Q] can be related to the phonon distri- 
bution in each medium if the heat flux q is expresed in 
terms of transmitted phonons : 

1 + 

q = — ^2v lx (p,k)huj(p,k)f 1 (p, k)t 12 (p,k) 

p,k 

+ y ^2 V2x(p, k)Hu)(p, k)f 2 (p, fe)*2i(p, k) (2) 

p , k 

where V is the volume of each medium supposed to be 
equal, Vi X is the group velocity in medium i projected 
along the direction x normal to the interface, fi is the 
mode dependant phonon distribution function in medium 
i, tij(k) is the wave vector dependant transmission coeffi- 
cient from medium i to medium j, and the sums runs over 
all polarizations indexed by p, and over wavevectors in 
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FIG. 2. (Color online) Steady phonon distribution function 
on both sides of the interface crossed by a heat flux. Red 
solid lines represent the local equilibrium distribution function 
faq{T(x)) given by the Bose-Einstein distribution eq. $Jfy. We 
have considered classical phonons for which the Bose-Einstein 
distribution is proportional to the temperature. Blue solid 
lines represent the out-of-equilibrium distribution / eq (T(a;)) + 
5f of incident phonons where Sf is given by eq. ©. We have 
considered a phonon mode propagating in the direction of the 
temperature gradient so that Sf < 0. We have also illustrated 
the graphical construction inherent to the definition of the 
equivalent equilibrium temperature (see the eqs. [13] and I14|l . 
Ai and A2 are the phonon mean free paths of the considered 
phonon mode in each medium. 



the first Brillouin zone corresponding to phonons cross- 
ing the interface, i. e those for which v\ x > and V2 X < 
respectively. In the following we will drop the variables 
(p, k) indexing the phonon polarization and wavevector 
to simplify the notations. We will refer to any quantity 
depending on (p, k) as mode-dependent. 
The problem of the determination of the interfacial con- 
ductance eq[T] relies on our knowledge of the phonon dis- 
tribution functions fi at both sides of the interface. The 
simplest reasoning is to assume that the phonons popu- 
lation can be described by the equilibrium distribution 
f eq given by the Bose-Einstein distribution: 



/eq(w,T) 



exp {huj /ksT) — 1 



(3) 



at the temperature Tj in the vicinity of the interface. The 
two sums appearing in eq.[2]can be contracted to a single 
sum on the phonon population coming from medium 1 if 
we invoke the principle of detailed balance in the situa- 
tion where the two media are at thermal equilibirum at 
the common temperature T\ when the flux q vanishes 11 . 
One arrives then at the Landauer formula for the inter- 
facial conductance 24 : 



G m = ^ hwv lx t 12 - ' '' 



1 

P j 



dT 



(4) 



Note of course that using the principle of detailed bal- 
ance, the Landauer conductance can be expressed as a 



function of the transport properties characterizing the 
medium 2: 



P . fe 



(5) 



The Landauer formula has commonly been used in the 
determination of the Kapitza conductance^^. Its limi- 
tations are well know n 11 ' 17 ' 19 ' 20 ' 27 : equation Q predicts 
a finite conductance when the two materials are identi- 
cal, i.e. when yk,t\2{k) = 1, which is of course contrary 
to the intuition: for an interface between similar media, 
the temperature drop should vanish whatever the flux 
g, leading to an infinite conductance 21 . Obviously, the 
problem is related to the use of two equilibrium distribu- 
tion functions in the flux cq(2] The previously mentioned 
paradox may be solved using the actual distribution func- 
tion consistent with the interfacial heat flow. This analy- 
sis has been done by Simons^ 7 - and generalized by CherJS 
and Landry and McGaughey 1 — . The out of equilibrium 
distribution function is supposed to obey the Boltzamnn 
transport equation (BTE) under the relaxation time ap- 
proximation 2 ^: 



dt 



fi fee 



(G) 



where / cq is the Bose-Einstein distribution given in eq. [3] 
and Tj is the mode dependant relaxation time supposed 
to depend only on the frequency U). In steady state, a 
solution of the BTE equation eq. [S] can be found in the 
form: 



fi(r f ) = f« l (T(f)) + 5f i {i f ) 



(7) 



where T(r) is the local value of the temperature. We as- 
sume in this way that the temperature is defined at any 
point of the material, an assumption which is reasonable 
if the phonon mean free path in each medium is not too 
large compared to the characteristic dimensions of the 
system, i.e. the distance between the interface and the 
heat reservoirs. Anyway, from a practical point of view 
in MD, one can always think of the local temperature as 
the mean kinetic energy of the atoms in a small volume 
encompassing the point r. If furthermore, we assume 
that in each medium the temperature profile is linear, 
an assumption which is again confirmed by NEMD sim- 
ulations^, then the deviation from the local equilibrium 
writes: 



Sfi(f) 



df e , 



Jcq 
-Ti-^-Vi 



dT 



VT 



(8) 



Hence, the excess of phonons propagating in each 
medium is proportional to the heat flux. Phonons trav- 
elling in the direction of the flux are in excess while 
phonons travelling in the opposite direction are depleted. 
A schematic representation of the distribution of incident 
phonons across the interface is displayed on fig. [5] Inject- 
ing the latter distribution function given by eqs. [7] and [51 
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in the interfacial flux q eq. [2 one arrives at: 



q = G oq (T 2 - Ti) - ^2 T ivl x hw^^t 12 ^\i 



p.k 



p.k 



where the temperature gradients are estimated on both 
sides of the interface. The two temperature gradients 
can be eliminated if we assume diffusive heat transport 



in each medium so that q — — Ai|^|i = — A2 1 2 where 



Xi denotes the thermal conductivity of medium i. 
interfacial conductance writes then: 



Gn 



1 - 013 - 02 



The 



(10) 



where we have introduced the fractions: 



012 = 



V 

p , k 



Tivl 



(11) 



and a similar equation for 02 1. The physical signification 
of 0i2 is clear : it is a measure of the fraction of the energy 
flux flowing across the interface that is transmitted. This 
coefficient varies typically between when all the phonon 
modes of medium 1 are reflected by the interface to 1/2 
when all the modes of medium 1 are transmitted. In 
particular, if we consider the case of similar materials, it 
is easy to show that the interfacial conductance eq. [10] 
diverges to infinity using the Peierls expression for the 
thermal conductivity^: 



fko 



dT 



012^12 - 1) 



(12) 



p . h 



where the last equality applies to the case of an interface 
which transmits all the phonon modes, V(p, k),ti2 = 1- 
Thus at least, eq. [TU] solves the paradox of the conduc- 
tance of the interface between identical materials. Note 
that we could have obtained the same expression for the 
conductance using the concept of equivalent equilibrium 
temperatures. By definition, the equivalent equilibrium 
temperature may be defined in a classical system in terms 



of the kinetic energy of the incident phonons in the vicin- 
ity of the interface. This condition is graphically illus- 
trated in fig. [5] and is mathematically expressed by : 



f eq (T^) = f cq (T 1 )+Sf 1 



(13) 



yielding 



Ti - Ai cos 9 



dT 

dx 



T(x — —Ax cos 9) (14) 



where the interface is supposed to be localized at x = 0. 
Here the 9 is the angle of incidence and Ai is the mean 
free path of the considered phonon mode. Thus, we have 
shown that the equivalent equilibrium temperature is the 
temperature of incident phonons at a distance of one 
mean free path away from the interface as exemplified 
in fig. [2j This explains why Aubry et al. obtained an 
expression similar to Simon conductance using the equi- 
librium distribution of phonons at a distance one mean 
free path from the interfac o 19 i 28 . Note that in the previ- 
ous discussion and in the formula used by Aubry et al., 
the equivalent temperature is a mode dependent quan- 
tity, as both 9 and Ai depend on the considered mode. 
Using the concept of equivalent temperatures may be 
dangerous because one may be tempted to believe that 
the phonon population is at equilibrium at a distance A 
away from the interface, which is of course wrong. It 
is nevertheless not surprising to find the same value of 
the interfacial conductance using the concept of equiva- 
lent temperature at a distance A, because the effective 
incident flux that may be transmitted by the interface 
comes from phonons which have not been scattered by 
other phonons before reaching the interface^ and as a 
first approximation if temperature gradients are not too 
large the corresponding phonon population may be de- 



scribed by /, 



cq V 



-Acosfl). 



In the following, it will be useful to express the different 
conductances in terms of the vibrational density of states 
(vDOS): 



5 P H = y^2 s ( 



£j r) 

p,k> 



(15) 



where the sum runs over the eigenmodes of the crystal in 
the first Brillouin zone. In the common case where the 
transmission coefficients depend only on the frequency uj 
and on the incident angle 9, the Landauer conductance 
is: 



tt/2 



ti2 (w,6>) cos 9 sin 9d9duj (16) 



and the fraction 12 becomes: 
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^ \ E p / " max 9l Au)Ti{u)\vi{u)\ 2 hw d J^ fj* tl2{uJ ,0) cos 2 flsinM?^ ( . _ ^ 

I E P Jo""' 1 5i , P )ri (w) I «i (u> ) | 2 fiw ^ 



where the factor 1/2 in the numerator comes from the 
integration over the azimuthal angle <fr and the integra- 
tion is carried out over the first Brillouin zone. oj max is 
the maximal frequency transmitted by the interface and 
its value will be discussed later and u)d,\ is the Debye 
frequency in medium 1. Again, a similar expression for 
the term /?2i can be obtained by permuting in the previ- 
ous equation the indexes 1 and 2. The challenge is now 
to specify the lifetimes Ti(uj) and the transmission coeffi- 
cients. We will discuss possible expressions for t\2 based 
on traditional interfacial transport models in the section 
IVII when we will analyze the conductance obtained by 
NEMD. 

So far, we have seen two formulae relating the interfacial 
thermal conductance to the energy transmission coeffi- 
cient: the Landauer fomula eq.|4] which assumes that the 
phonons on both sides of the interface are at equilibrium, 
and the general formula eq.fTOl which accounts for the ac- 
tual out of equilibrium distribution of the phonons in the 
vicinity of the interface. Now the question that we want 
to answer is: what do we measure in a molecular dynam- 
ics simulation? Intuitively, in NEMD simulations where 
the system is crossed by a flux, we should measure a con- 
ductance given by eq.lTOl because the system is subject to 
a large temperature gradient (on the order of 1 K/nm !) 
and the phonons can not be considered locally at equilib- 
rium. On the other hand, it seems reasonable to consider 
that in an equilibrium simulations where thermally in- 
duced fluctuations of the interfacial flux are probed, one 
should measure the Landauer conductance eq. 0] rather 
than the non-equilibrium conductance eq. 1101 We will 
make this point more quantitative in the next section. 



III. GREEN-KUBO FORMULAE: 
CONDUCTANCE FROM EQUILIBRIUM 
FLUCTUATIONS 

In this section, we derive Green-Kubo formulae for the 
interfacial conductance. We will prove that the Puech 
formula traditionally used in equilibrium simulations to 
measure the interfacial conductance is exactly given by 
the Landauer conductance eq. (|4]) and thus differs from 
the non-equilibrium conductance G neq (eq. (fTU)0 . We will 
also propose an equivalent formula easier to evaluate in 
molecular simulations. 

The general idea behind Green-Kubo formulae is that 
the regression of the fluctuations of an internal variable 
obeys macroscopic laws. In the case of interfacial heat 
transfer, the relevant variable is the interfacial flux q and 



the corresponding Green-Kubo formula reads: 

Gpucch = (q(t)q(0))dt (18) 

This formula has been used for solid/liquid interfaces lD 
and superlattices as wel l 16 i 29 . However practically in a 
MD simulation, the expression of the heat flux q involves 
only atoms near the interface^ and it is also sometimes 
difficult to estimate the plateau in the heat flux corre- 
lation function in eq. 1181 In the following, we will show 
that we can improve the statistics on the determination 
of the interfacial conductance by measuring the fluctua- 
tion of the mechanical energy of each solid on both sides 
of the interface. In passing, we will show that for the 
case of solid/solid interfaces, the Puech formula eq. (fT8|) 
identifies with the Landauer conductance eq. Q. 
To this end, we consider two semi-infinite media sepa- 
rated by an interface whose area is denoted A. The 
two media are supposed to be at thermal equilibrium 
at the same temperature T, and the energy in each 
medium can change only because of exchange of energy 
with the other medium through the interface. Generally 
speaking, the energy fluctuation in each medium is^S: 
(SEf) = k B T 2 /(l/C v i + l/C v2 ) where E t (t) is the instan- 
taneous mechanical energy of the medium i, and C v \ , C V 2 
are the specific heat characterizing the two media. Note 
that the relevant statistical ensemble to describe the fluc- 
tuations of Ei is neither NVE because only the total 
energy E\ + E% is conserved nor NVT because stricly 
speaking each system is not in contact with a thermostat 
but with a system which is comparable in size. We refer 
the reader to Stephenson22 for a derivation of the fluc- 
tuations of the different quantities in this situation. The 
classical NVT formula is however recovered when one of 
the two media (say 2) has a large number of degrees of 
freedom so that C V 2 ^> C v ±. In the following, we will 
assume that the two media have the same specific heat 
so that the energy flucuation in each medium is: 

(SEf) = k B T 2 C v /2 (19) 

This hypothesis will not affect the final result but al- 
lows to simplify the notations all along the derivation. 
The fluctuations of the interfacial flux q are related to 
the fluctuations of the energy in medium 1 through the 
energy conservation equation: 

where q is the instantaneous value of the interfacial en- 
ergy flux flowing from the medium 1 towards medium 2, 
which in the situation considered fluctuates around zero. 
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This flux may be expressed in terms of excess phonon 
occupation number Srif. 

I + 

g = y / ] Vx x tvj)8nitx2 + } J viJ\us&nitz\ (21) 



where the excess phonon occupation number is simply re- 
lated to the fluctuation of the energy: SEi = fuxibn i g 

and here k is a shorthand notation representing the 
wavevector and the polarization. In the following, it will 
be useful to rewrite the energy conservation: 



~~dT 



1 



frwdrii 



1 v ■> hujSn 2 



k k 



(22) 



where we have introduced the mode dependent relaxation 
times: 

= V/ Avi x t 12 if V\ x > (first sum in the rhs of eq. I2"12j) 

= V^/^|w 2 a;|*2i if < (second sum) (23) 

These relaxation times may be interpreted as interfacial 
scattering terms and are independent of the bulk phonon 
relaxation times. If we assume the different modes to be 
independent and consistently with eq. (|19[) characterized 
by a variance 



k B T 2 c v 



(24) 



where c v — fiw^^pp is the mode dependent specific heat, 
the energy autocorrelation function follows : 



(SE^SE^O)) = 



fc 



— exp(-|t|/r £ ) 



(25) 



where we have used the total energy conservation 8E\ = 
—SE2- Differentiating this latter equation, one arrives at: 

,dEi(t) . , v-^ k B T 2 c v . . . . . , . 

(—7^^1(0)) = - E -^T^W ex P (-1*1/^) 



dt 



2 ^ 



k 



k B T 2 c v 5{t) 



(26) 



where sgn(i) is the sign function and the second term in 
the right hand side comes from the discontinuity of the 
derivative of exp(— |i|/rg) at the original The sums over 
all the wavevectors k may be transformed in a sum run- 
ning over the modes crossing the interface if we express 
the detailed balance condition: 



C v Vl x tl2\v 1 x>0 

yielding for t > 0: 



~C v V2xt21 \v 2 x<Q 



(27) 



{ ^Ml 6El (0)) = -i]T^fc s T 2 c„i 12 « 1:c exp(-i/^) 

fc 

+ 

+ 2j2 k BT 2 c v 8(t) (28) 



Note that in the thermodynamic limit the first term van- 
ishes. For a finite system and when the time t — > + , 
the second term involving a Dirac distribution may be 
neglected and one has : 



1 



f dC EE \ 



Ak B T 2 \ dt J t=Q i 



I + 

— 22 c v t 12 v lx = G cq (29) 



where Cee^) = {8E\{t)8E\{Q)) is the energy auto- 
correlation characterizing solid 1 and G eq is the Lan- 
daucr conductance defined in eq. [4] The previous equa- 
tion eq. (|29j) is a new Green-Kubo formula for the in- 
terfacial conductance which relates the slope of the en- 
ergy autocorrelation function at the origin to the Lan- 
dauer conductance. We will show in the next section 
that this formula may be easier to estimate in a MD 
simulation than the classical Puech formula which re- 
quires to identify a plateau in the running integral of 
the heat flux autocorrelation function. Alternately, we 
can also relate the Puech formula to the Landauer con- 
ductance in the thermodynamic limit by remarking that 

.d5n 1 ^(t) d5n 1 g(0) 



dt 



dt 



) = -^(S ni M)Sn ui (0)) to 



arrive at 



,dEx(t) dEi(0) v 



dt 



dt 



Ek B T 2 c v 
— - — expi 



1 



l*l/t) 



2 MBT^Y j c v t 1 2V lx 6{t) (30) 



In the thermodynamic limit, in principle the first term 
cx A(A/V)22- is negligible compared with the second oc A, 
and one has the following Green-Kubo equation: 



1 



Ak B T 2 



,dE(t) dE(0) 
"~dt dl~ 



dt = G, 



cq 



(31) 



This equation is exactly the Puech formula used to cal- 
culate the liquid/solid Kapitza resistanc e 15 ' 38 . We have 
shown that for solids/solids this formula identifies with 
the Landauer conductance. It is important to realize 
that the previous formula has been derived for an infi- 
nite system size. For a finite system on the other hand, 
the running integral 

' ,dE(t') dE(0) , 1 



Ak B T 2 J N dt dt 



dt' = — c v t 12 vix 



I + 

-y22, Cvtl2Vlx ( 1 ~ ex v{-t/T%)) (32) 



will consist of two parts: the first term is the Lan- 
dauer conductance, the second term is negative and cor- 
responds to the final decay of the running integral. 
In the next section, we will compare the formulae 
eqs. Q and (JTSJ) to the results of NEMD simulations. 
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IV. SIMULATIONS 
A. Lennard-Jones systems 

We now study how the previous formulae may be 
used in molecular dynamics simulations to estimate the 
Kapitza conductance between two solids. All the follow- 
ing results have been obtained for the case of the inter- 
face between Lennard-Jones solids. There are numerous 
advantages to work with LJ solids. The first is the sim- 
plicity of the interaction potential as compared to many 
body potentials used to model semi-conductors. This has 
an important practical consequence as it allows to run 
simulations with large system lengths because of the rel- 
atively short computational times required. Also, from 
a thermal point of view there is no need to worry about 
optical phonon modes. 



B. Structures 

We will consider systems consisting of two perfect fee 
Lennard-Jones solids whose interface is orientated along 
the crystallographic [100] direction. The section of the 
system is fixed to 6ao X 6ao where ao is the fee lattice 
constant, and the thickness of each medium has been 
varied between 10 and 50 ao. A typical initial con- 
figuration is represented in figure [3] All the atoms of 
the system interact through a Lennard-Jones potential 
= 4e((tr/r) 12 — (er/r) 6 ) truncated at a distance 
2.5a. A single set of energy e and diameter a charac- 
terizes the interatomic interaction potential. As a re- 
sult, the two solids have the same lattice constant ao, 
and the interface may be considered perfect. To intro- 
duce an acoustic mismatch between the two solids, we 
have considered a mass mismatch between the masses 
of the atoms of the two solids, characterized by the 
mass ratio m r = m^/mx, which will take typical val- 
ues between 1 and 10. From now on, we will use real 
units where e = 1.67 10~ 21 J; a = 3.4 10~ 10 m and 
mi = 6.63 10~ 26 kg, where these different values have 
been chosen to represent solid Argon. With this choice 
of units, the unit of time is r = ^Jma 2 Je = 2.14 ps, 
the unit of thermal conductivity is A = ^yzjrn ~ 18.8 
10~ 3 W/K/m and the unit of intcrfacial conductance is 
G = k B /(Ta 2 ) ~ h&MW/K/m 2 . The different interfaces 
have been prepared as follows: first the structures have 
been generated by mapping the space with fee structures 
using the lattice parameter of the fee LJ solid at zero 
temperature 3 ^: ao(T = OK) = 1.5496c. The structures 
have been then equilibrated at the final finite tempera- 
ture T = 40 K using first a Berendsen thermostat and a 
barostat at atmi£. Once the instantaneous temperature 
has increased to a value close to the final expected tem- 
perature, we have switched off the Berendsen thermostat 
and used a Nose Hoover thermostat. The total equilibra- 
tion time lasts one million time steps which correspond 




FIG. 3. (Color online) Configuration studied: interface sep- 
arating two Lennard Jones fee solids having the same lattice 
constant but different masses. 



to a total time of 4, 28 ns. All the systems studied have 
been equilibrated at the temperature T = 40 K, and the 
lattice parameter at this temperature has been found to 
be: ao = 1.579(7. In EMD, periodic boundary conditions 
have been applied in all spatial directions so that the sys- 
tem represented is a superlattice 2 -. On the other hand 
in NEMD we use periodic boundary conditions only in 
the directions parallel to the interface. 



Computing the interfacial conductance with 
NEMD 



Alternately, we will compare the results of EMD to 
NEMD. The principle of these latter simulations has been 
already described elsewher e) 11 ' 34 and we just focus here 
on the details of the technique. We impose a thermal 
flux perpendicular to the interface between the two solids 
by thermostatting in each medium two layers of atoms 
remote from the central interface at the respective tem- 
perature T c = 40 - 3.6K and T H = 40 + 3.6 K, while the 
end atoms are maintained fixed. The size of the cold and 
hot regions has been found to have negligible effect on 
the measured conductance. After a number of time steps 
varying between 500000 for the smallest system to 5 mil- 
lion for the largest, we monitor the temperature profile in 
each medium using the kinetic energy of the particles. A 
typical example of the corresponding temperature pro- 
file is shown in figure |4] zooming on the vicinity of the 
interface. The interfacial conductance is obtained from 
the heat flux and the temperature jump across the in- 
terface, these latter quantities being measured using the 
heat power delivered by the heat source which is moni- 
tored during several million of time steps. The tempera- 
ture jump AT is obtained by extrapolation of the linear 
profiles in the two media as shown in figure 01 In the 
following, we will present results for the interfacial con- 
ductance obtained using 5 independent simulations. 
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FIG. 4. (Color online) Stationary temperature profile across 
the interface obtained by NEMD. The symbols are the sim- 
ulation data while the solid lines are the linear extrapolation 
used to compute the interfacial temperature drop AT. The 
position of the interface is located by the vertical dashed lines. 



D. Computing the interfacial conductance with 
EMD 



As we have discussed in the previous section, there are 
several formulae to compute the interfacial conductance 
from EMD simulations probing the energy flux between 
the two solids. First, we will consider the energy auto- 
correlation formula in eq. To obtain the value of the 
Landauer conductance, one needs to compute the time 
derivative of the corresponding energy autocorrelation 
function. To this end, we have recorded the instanta- 
neous value of the mechanical energy E™ of each solid: 



rav] 
2 J 



E v(f 

j,kei 



(33) 



where the first term represents the total kinetic energy of 
the solid i and the second is the potential energy between 
atoms belonging to the solid i. Note that this definition is 
somewhat a little bit arbitrary and we could have chosen 
to include in the mechanical energy the cross interaction 
term J2jetkej ^(^f ~ ^0- We have not observed signif- 
icant differences in the value of the Landauer conduc- 
tance as compared to the first definition. To determine 
the value of G eq , we have computed the energy autocor- 
relation function (EACF) in each medium. The instan- 
taneous value of the mechanical energy of each solid has 
been recorded every two time steps in the course of long 
NVE simulations corresponding to a total of 1 million 
time steps. The EACFs have been obtained by averaging 
over 10 initial independent configurations. An example 
of the averaged EACFs is shown in fig. [5] The EACFs rel- 
ative to the two solids are practically indistinguishable. 
Note the existence of very small oscillations of the EACFs 
at long correlation times. Given the value of the mean 
sound velocity in the system c ~ 1.2 nm.ps -1 , these os- 
cillations should probably correspond to long wavelength 



phonons which have travelled ballistically across the sys- 
tem several times, thus creating "echoes" in the correla- 
tion functions. To estimate the value of the time deriva- 
tive at time t = + appearing in eq. 1291 we have fit the 
EACFs with a single exponential function between a time 
t ~ 5 ps and up to a time where the EACF has decreased 
by a factor 10 as compared to the initial value. Using 
the fit Cee^) = Cee(0) exp(— t/r), the conductance is 
G = 2^f^2 where C EE (0) is found to be k B T 2 Vpc v /4 
to a good approximation and the factor 2 in G comes 
from the fact that there are two interfaces due to the pe- 
riodic boundary conditions. The uncertainty in the de- 
termination of the value of G is found to be typically 20 
percent for 10 independent configurations and of course 
it decreases with the number of realizations of the sys- 
tem. Finally, we want to emphasize that we have ob- 
served that for large systems, the EACF decreases very 
slowly in good agreement with the previous mode analy- 
sis eq. (|23p which predicts that the mode relaxation times 
scale as the system length. 

Alternatively, we have also analyzed the conductance us- 
ing the Puech formula eq. ([IB")) where the instantaneous 
value of the flux may be estimated in the course of a MD 
simulation using the power of the interfacial forces^: 



q = 



E 



' Fir. 



(34) 



Note that this expression of the flux q involves only atoms 
in the vicinity of the interface, while all the atoms of the 
system contribute to the expression based on eqs. (j2T))) 
and (|33[) . Figure [6] displays the running integral in the 
Puech formula eq. ([TBI) calculated using simulations for 
the same system considered in figure [5j The two curves 
correspond to the two interfaces of the system (remember 
the periodic boundary conditions) . The running integrals 
display first a peak and then slowly decrease. Note the 
echoes in the upper curve. We have found that it was 
difficult to define unambiguously a plateau eventhough 
we have considered here an average over 30 independant 
configurations. The slow decrease has been also observed 
in the determination of liquid / solid conductance 1 - and is 
predicted in eq. (|32[) . Indeed it is a common problem for 
a finite ergodic system that the Green Kubo formula pre- 
dicts a vanishing transport coefficient^! and in practice 
the running integral should be estimated at an interme- 
diate time To where the integral has not yet significantly 
decreased. The problem in heat transfer simulations of 
solids is that the spectrum of relaxation times spans 
several decades and defining an intermediate time To is 
not obvious in this situation. This difficulty is somewhat 
circumvent in the formula eq. ([29} as it does not require 
to estimate a plateau. 

Finally, we compare the value of the interfacial con- 
ductance to the expression proposed by Rajabpour and 
Volz 14 for a classical system : 



1 

G 
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+°° (ST(t)ST(0)) f 11- 
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(35) 
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FIG. 5. (Color online) Energy autocorrelation functions of the 
two fee Lennard- Jones crystals separated by a planar inter- 
face obtained with molecular dynamics simulations. Dashed 
lines show the exponential fit. The parameters are: Total 
Length=40 ao; T = 40 K; mass ratio mr = 2. 

where N\ and N2 are the number of degrees of free- 
dom characterizing each medium. Practically, the inter- 
facial resistance is obtained by fitting the kinetic energy 
autocorrelation function with a single exponential hav- 
ing a decay time r. The conductance is then given by 
G = AkBpV/{2r). Figure [7J displays the kinetic energy 
autocorrelation function obtained by averaging over 10 
independent simulations for the same system as consid- 
ered before. As noted before^, the kinetic energy dis- 
plays a first fast decrease followed by a longer decrease, 
which is fitted with a single exponential with a relaxation 
time r. This latter time is used to obtain the interfacial 
conductance G — AkupV j '(2t). Note the oscillations in 
figure [7] due to the conversion between kinetic and po- 
tential energy. These oscillations occur with the same 
period than the period of echoes observed in the energy 
correlation function fig [5] The value of the interfacial 
conductance obtained G — 47 ± 12 MW/K/m 2 is smaller 
than the value obtained with the energy correlation func- 
tion G = 60 ± 12 MW/K/m 2 but within the error bars. 
Hence, the two methods give consistent results and com- 
parable error bars. 

We conclude by saying that compared to the Puech for- 
mula eq. (fT5)) , the new Green-Kubo formula eq.[55]is eas- 
ier to evaluate in a MD simulation because: 1. we do not 
need to estimate a plateau in a running integral; 2. the 
new formula involves all the atoms of the system while 
the Puech formula involves only atoms in the vicinity of 
the interface. As a result, the statistics is improved. 



V. FINITE SIZE EFFECTS 

In this section, we compare the finite size effects in the 
determination of the conductance using both EMD and 
NEMD. Figure [5] a. displays the length dependance of 
the Kapitza conductance obtained with equilibrium sim- 
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FIG. 6. (Color online) Interfacial conductance of the same 
interface as fig. [5]calculated using the Puech Green-Kubo for- 
mula eq. 1181 The two curves correspond to the two interfaces 
of the system. Each curve is an average over 30 independant 
trajectories. Same parameters as figure [5] 
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FIG. 7. (Color online) Kinetic energy autocorrelation func- 
tion for the same system as fig. [5] 



ulations Gemd- It is found that Gemd decreases with the 
system length. We have not studied the conductance of 
systems longer than 100 ao because as explained above it 
leads to very long relaxation times (see eq |2"5)) and the 
determination of the equilibrium conductance becomes 
costly. In figure [8] b., we quantify the finite size effects on 
the conductance obtained with NEMD. The values ob- 
tained are consistent with the results of Stevens et alJ^. 
Note the values of the NEMD conductances which are 
larger than the EMD conductance by a factor 5! This 
discrepancy will be analyzed in detail in the next sec- 
tion. We focus our attention here on the less severe size 
effects displayed by the NEMD conductance as compared 
with the EMD. The finite size effects in the EMD m ethod 
are quantitatively analyzed in the appendix IVIIII 

The general idea is the following: in the EMD simu- 
lations, there are two interfaces between the two media 
to be considered because of the periodic boundary con- 
ditions as sketched in fig. [SJ These two interfaces are not 
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FIG. 8. (Color online) Interfacial conductance obtained by 
EMD (top figure) and NEMD (bottom figure) as a function 
of the length of the system in units of the crystal monolayer 
ao ■ The solid red line on the top figure displays the theoretical 
formula eq. (|36[) . Same other parameters as fig. [5] 



necessarily independent from a thermal point of view: 
longwavelength phonons having a long mean free paths 
can create correlations between the instananeous value 
of the thermally induced flux at two adjacent interfaces. 
More precisely, if we denote by A and A' the two inter- 
faces, the energy conservation writes: = qA + qA' and 
the calculation of the equilibrium conductance involves 
cross terms of the form (<?a(*)<w(0)) and {qA> (t)qA{0)), 
while the "intrinsic" interfacial conductance is given by 
the term: (q A (t)qA(0)) = (q A >{t)qA>{0)). Clearly, the 
cross terms will be relatively important at small interfa- 
cial separation L/2 because a majority of phonons modes 
will have a mean free path larger than L, while they 
should vanish in the limit L — > oo. These cross terms are 
quantified in the appendix, under the assumptions of the 
interface between Debye solids with a constant transmis- 
sion coefficient t\2 1 an assumption assessed a posteriori 
as shown in the next section fVI El where we will show that 
the EMD results are well described by the DMM model. 
We have also assumed that the phonon relaxation times 
are described by the Callaway model that we will discuss 
in the next section (cf eq. I50[) . The prediction derived in 



the appendix I Villi may be written: 



G{L) = Goo ( 1 + c | - 



3/2N 



(36) 



where Goo is the conductance caracterizing the interface 
between semi-infinite media, c is a numerical constant 
and 



h 
G, 



(37) 



is the phonon correlation length in medium i where 
Gu = ^riikBCi is the EMD conductance between two 
identical media having the properties of medium i (see 
also the next section). In principle, one should define two 
phonon correlation lengths characterizing the two media, 
but in the case of mass-mismatch Lennard- Jones solids, 
the correlation length is the same for the two media and 
is £ ~ 6.3ao at T — 40 K where we have used the value of 
the thermal conductivity obtained by Green-Kubo simu- 
lations by McGaughey and Kaviany^. Figure [8] a. com- 
pares the EMD data to the theoretical expression eql3"rJl 
where we have obtained that the constant c ~ 5 above 
the predicted value c = 2.75 (see the appendix lVIII[) . The 
disagrement may be due to the use of the DMM model 
which as we will see in the next section tends to over- 
estimate the conductance obtained in EMD thus under- 
estimating the constant c. Note however that we could 
have fit with the same accuracy the EMD data using a 
functional form G(L) = Goo(l + C (£,/L)) and that the 
G(L) ~ L~ z / 2 scaling comes in our analysis from the 
Callaway assumption. What is important to remember 
is that the EMD conductance decays algebraically with 
the system length with a characteristic length propor- 
tional to the phonon correlation length £. On the other 
hand in NEMD the length dependence is smaller because 
the distribution of phonon mean free paths is cut due to 
the presence of heat resevoirs. 



VI. COMPARISON BETWEEN THE EMD AND 
NEMD CONDUCTANCES AND THEORETICAL 
MODELS 

In this section, we will compare the values of the inter- 
facial conductance obtained either by EMD and NEMD 
(and extrapolated to infinite system lengths) to the ex- 
pression of the conductance eqs. (dJ) and (fit))) where we 
should specify the value of the phonon transmission co- 
efficient ti2- To this end, we will consider two classical 
models for interfacial phonon scattering: the AMM and 
the DMM. We will generalize these two models to de- 
scribe the non-equilibrium conductance cq. (|10p . We will 
present in passing usefull approximate analytical expres- 
sion to estimate both the Landauer conductance and the 
general conductance combined with the AMM model. 
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FIG. 9. (Color online) Postulated origin of the length de- 
pendence of the conductance Gemd measured in EMD: long 
wavelength phonons may travel ballistically between the two 
interfaces A and A' thus creating thermal cross-correlations 
as measured by (qA{t)qA'(0)}- This figure displays also the 
notations used in the appendix I Villi to quantify this effect. 



A. Debye approximation 

All along this section,we will make the assumption of 
Debye solids. In the Debye approximation, the solids 
are assumed to have a constant group velocity which 
depends on the polarization mode — . Most often an 
additional assumption is made consisting in assuming 
the same acoustic velocity for each polarization^. For 
a three-dimensional crystal, this latter is defined as : 



Ccff 



2c T + c L 



(38) 



where the indexes T and L refer to the transverse and 
longitudinal polarizations respectively. Under this as- 
sumption, the vDOS is : 



2tt 2 c 



(39) 



In the following we will drop the subscript " efi" and char- 
acterize the averaged sound velocity in medium i by Cj. 



B. Acoustic Mismatch Model 

In the AMM, the phonons traveling towards the inter- 
face are assumed to see the interface as a sharp disconti- 
nuity of acoustic impedance Z% where 



which strictly speaking holds as long as the incident angle 
is smaller than the critical angle 9 C — arcsin(c2/ci). Here 
0\ and 62 are the incident and refraction angles respec- 
tively. Above the critical angle, as for the electromag- 
netic waves, internal reflection occurs and the incident 
phonons are totally reflected. For the Si/Ge interface 
for which the ratio of the acoustic velocities is approx- 
imately 1.5, the critical angle is 6 C ~ 40deg and a sig- 
nificant fraction of phonons are totally reflected by the 
interface. We have also assumed that a phonon conserves 
its polarization, i.e. there is no mode conversion and c\ 
and C2 denote the acoustic velocities in media 1 and 2 
corresponding to the same polarization state. Another 
assumption behind eq. [41] is that the scattering is elas- 
tic, i.e. refracted and reflected phonons conserve their 
frequency. As a consequence, phonons having frequency 
above the Debye frequency of the softer solid are confined 
in the hard solid and not transmitted by the interface, 
i.e. the transmission coefficient is supposed to vanish. 
For phonons having frequencies smaller than the debye 
frequency of the softer solid, the transmission coefficient 
is derived from the Snell law2£ eg ATI 



t 12 (0J,fXl) 



4ZiZ 2 /ii^2 
(Zxfii + Z2II2) 2 
otherwise; 



uj < min(cjDi, ^02) 



(42) 



where we have introduced jii = cos (9.;, and again it is 
implied that the incident angle is smaller than the crit- 
ical angle. At high temperatures, the regime relevant 
to classical molecular simulations where the equilibrium 
Bose-Einstein distribution / eq (w) — > ksT '/Huj, the AMM 
conductance which is calculated using the Landauer ex- 
pression eq. |4] may be written : 



G- 



AMM 



-ni&BCi 



^12(^1)^1^1 



(43) 



(40) 



where n\ denotes the number density of medium 1. We 
have supposed without loss of generality that the medium 
denoted 2 has the lowest Debye frequency. The factor 
(— ) comes from the phonon confinement of high fre- 
quency phonons in medium 1. The AMM conductance 
eq. (|4"3"f should be evaluated numerically. Alternatively, 
one can obtain tractable analytical expressions for the 
AMM conductance, if we assume that when the acoustic 
contrast between the two solids is large, the transmis- 
sion coefficient t\2 is dominated by phonons propagating 
with a small refraction angle, i.e. ^2 — !• Under this 
approximation, the AMM conductance is given by the 
approximate form: 



is given by the product of the mass density pi by the 
acoustic velocity in the medium i. As a result, phonons 
may be reflected by the interface or refracted on the other 
side of the interface following the equivalent of Snell laws : 



sin c/i sin t/ 2 



ci 



C-2 



(41) 
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-nxk B Ci 



j-appx 



(44) 



where J^ ppx depends on the acoustic ratio f3 = Z2/Z1: 



1 i 
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(45) 
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As shown in the Appendix IIX1 the approximation eq (|44"]) 
gives a very good description of the AMM conductance 
over a wide range of acoustic contrast. 



C. Diffuse Mismatch Model 

The previously described AMM model is supposed to 
predict the transmission of phonons of large wavelengths 
which behave as plane waves experiencing specular re- 
flection or refraction at the interface. This model is 
commonly thought to apply at low temperatures where 
only long wavelength phonons are populated. At higher 
temperatures, interfacial scattering is thought to be dif- 
fuse like essentially because a majority of phonons have 
wavelengths comparable or even smaller than the inter- 
facial roughness. This idea motivated the development 
of the DMM introduced by Swartz and Pohl2J£ which 
assumes that the phonons experiencing scattering at the 
interface loose totally the information about the medium 
where they come from. As a result, the probability that 
a phonon experiences a reflection in medium 2 is equal 
to the probability that a phonon is transmitted from 
medium 1 towards 2 : 



where A4 is a material parameter which depends on 
the temperature. Under this assumption and if inter- 
facial scattering is supposed to be specular, the non- 
equilibrium conductance takes the form : 



G 



<-T„, 



1 _ f (ff") (lo ^1*12(^1)^1 + tj- Jo AtiM2*i2(Ati)d/ii 

(51) 

w here /.12 denotes the c osine of the refracted angle 37 :/i2 = 
y/l — (c2/ci) 2 (l — /if). Again, the conductance eq. (j5"Tj) 
can be approximated: 
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(52) 



where 



jappx j g j e g ne( j m e q.(g5]) anc i : 



/3 3 



+ 3/3 log ( — — 



(53) 

The accuracy of the approximation eq (|52[) and a hner 
approximation are presented in the Appendix IIXI The 
conductance obtained using the DMM transmission co- 
efficient is : 



tl2 — 1 — *21 



(46) 



for the particular mode considered. Writing the total flux 
in medium 2 together with the previous amnesia condi- 
tion yields the transmission coefficient : 



h2(u) = 



C2ff2M 

ciffi(w) + 



(47) 



and as for the AMM model, it is implicitely assumed 
that high frequency phonons are confined in the harder 
material: 



G 



DMM 



-nik B 



4 + (ci - C2 ) 2 



(54) 



Again we note that when the two media are similar, 
ci = C2 and the previous equation for the conductance 
predicts a finite conductance G(ci = C2) = |ni/csci. 
This new paradox can be traced back to the use of the 
DMM transmission coefficient eq. (j47|) which tends to- 
wards 1/2 when c\ — > c\. This problem disappears using 
the AMM transmission coefficient because the denomi- 
nator of eq. (|5"Tj) tends towards when the two media 
are identical. 



£12(^0=0 if w > min(u>D, 1,^73,2) (48) 

Since the transmission coefficient doesnot depend on the 
incident angle, the DMM conductance has a simple ex- 
pression: 

/iDMM 3 C 3 , . 

G eq = -nifcs-2— 2 (49) 
t 4T4 

D. Generalized conductances 

To obtain tractable expressions for the non-equilibrium 
conductance eq. ()10j) which depends on the fractions fiij 
eq. (|11|) . we need to do an hypothesis regarding the 
frequency-dependence of the phonon lifetime Ti{ui). The 
simplest is to assume that the phonon lifetime 7$ is con- 
trolled by Umklapp processes obeying Callaway model 22 : 

Ti(u) = A,lj- 2 (50) 



E. Interfacial conductance of a series of 
mass-mismatched Lennard-Jones solids 

In this subsection, we compare the conductances ob- 
tained using both EMD and NEMD simulations to the 
previous equations for the interfacial conductance, re- 
spectively given by the AMM model eqs. (|i5|). the DMM 
model eq. (|4"Tj|) and the generalizations eqs. (|5"Tj) . and 
eq. (|54l) . In evaluating these different expressions for 
the case of the interface between Lennard-Jones solids, 
we have used the values of Argon: c\ = 1250 m.s -1 for 
the average sound velocity of the harder medium and a 
number density n = 2.57 1 28 m _1 . 
In figure 1 101 we have reported the values obtained us- 
ing EMD and NEMD simulations for the interfacial con- 
ductance characterizing the interface between LJ solids 
having a variable mass ratio. This ratio has been varied 
between 1 and 10 so as to change the acoustic impedance 
ratio Z\ /Z2 — \pn\Jni2 between the two media between 
1 and 0.3. The EMD values have been obtained using 
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the finite size scaling analysis described before and the 
extrapolation to infinite system length as described in 
the previous section [Vj On the other hand the NEMD 
values have been obtained using a total system length 
of 200 a . The trend displayed by the NEMD data is 
very similar to the NEMD simulation results of Landry 
and McGaughey for the Si/heavy Si interface^. Strik- 
ingly and as already in the previous section, the EMD 
and NEMD values may differ significantly depending on 
the acoustic contrast between the two solids. In par- 
ticular, when the dissimilarity between the two solids is 
small, the NEMD conductance is larger than the EMD 
value by more than one order of magnitude ! Note that 
the corresponding impedance ratio ~ 1.5 are typical of 
AlAs/GaAs interface? 3 . Even for dissimilar solids like 
Si/Ge for which the impedance ratio is ~ 1.7, the dif- 
ference may reach a factor 3! This discrepancy may be 
simply explained: as we showed, the EMD conductance 
yields the Landauer expression of the conductance eq. Q 
while the NEMD value should be akin to the Simon con- 
ductance eq. (|10p . The difference between the two values 
of the conductance is quantified by the fractions of "out- 
of-equilibrium" phonons j3±2 and P21 (eqs. lll[) which tend 
to make the denominator of eq. (|10[) vanishing when the 
acoustic properties of the two solids become comparable. 
In this limit, the difference between the general expres- 
sion eq. (|10p and the Landauer conductance may be very 
large, yielding the divergence of the NEMD conductance 
when the two solids are similar. In figure I1Q[ we have 
also compared the EMD values to the AMM and DMM 
models which are consistent with Landauer formalism. 
Based on the analysis of the conductance at the inter- 
face between similar solids, we conclude that the DMM 
model gives a relatively good description of the EMD con- 
ductance, while the AMM model overpredicts the EMD 
values by a factor 2. Note however that the difference 
between the AMM and DMM models is not that large 
for dissimilar materials. The small discrepancy between 
the simulation values and the DMM model may come 
from our assumption of Debye solids in a situation where 
a fine description of high frequency modes is required, as 
the DOS of the two solids strongly overlap and the max- 
imal frequency transmitted by the interface w max tends 
towards the Debye frequency of the harder solid. Regard- 
ing the NEMD values, it is clear that the generalization 
eq. (|5ip based on the acoustic transmission coefficient de- 
scribes quite satisfactorily the divergence of the NEMD 
conductance. Equation (pyT| which relies on a diffusive 
transmission coefficient underpredicts the NEMD con- 
ductance by a factor larger than 5 for typical values of the 
acoustic impedance ratio. This is not completely surpris- 
ing since as we discussed before, if interfacial scattering is 
diffuse, the interfacial conductance does not diverge when 
the two solids become similar. Also importantly, we have 
seen that interfacial phonon transmission in EMD simu- 
lations is controlled by diffuse events, while it becomes 
determined by the acoustic properties of the two solids 
when a thermal flux is imposed. Hence, we conclude 
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FIG. 10. (Color online) Interfacial conductance determined 
by EMD and NEMD as a function of the mass mismatch be- 
tween the two solids. The simulations results are compared 
with the different theoretical expression AMM (eqESJ, DMM 
fca I49[) and the generalizations eqs. (|51[) and (|54[1 . The tem- 
perature is T = 40 K. 



that the energy transmission coefficient is not an intrin- 
sic property of an interface, and it may depend on the 
nature of the source of thermal flux (i.e. external heat 
reservoirs vs. internal fluctuations). Given the results 
of the simulations, we are tempted to conclude that in 
equilibrium simulations, thermal fluctuations destroy the 
correlations between incident and transmitted phonons 
so that the amnesia condition eq. (j46|) is verified and the 
conductance is well predicted by the DMM model. In 
particular when the two media are similar, one recovers 
the fact that a phonon in excess will have a probability 
1/2 to be transmitted and 1/2 to be reflected, which is 
consistent with the EMD values obtained in this limit. 
On the other hand, in a NEMD simulation the situation 
is quite different: indeed phonons travelling across the 
interface see the interface as a sharp discontinuity which 
creates strong correlations between incident and trans- 
mitted phonons. Because the thickness of the interface is 
smaller than the phonon wavelengths, the transmission 
and reflection coefficients will be in these conditions con- 
trolled by the acoustic impedances of the two media, and 
in the limit of similar solids the transmission coefficient 
should approach 1. This may explain the difference in 
transmission coefficients between EMD and NEMD sim- 
ulations. 



VII. CONCLUSION 

In conclusion, we have analyzed two methods to 
measure the thermal Kapitza conductance between di- 
electrics using molecular dynamics. We have proposed a 
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new Green-Kubo formula (eq. [29l) to measure the inter- 
facial conductance using equilibrium EMD simulations. 
This formula is easier to evaluate in a molecular dynamics 
simulations as compared to the classical formula eq. (|T5]) 
because it avoids to estimate a plateau in the running 
integral of a correlation function. Also the statistics is 
improved because the new formula involves all the atoms 
of the system when the Puech formula considers only the 
atoms in the vicinity of the interface. We have also an- 
alyzed finite size effects in EMD and showed that their 
origin is the correlation between the interfaces created 
by long wavelength phonons which travel ballistically 
across the periodic simulation cell. On the other hand, 
in NEMD the distribution of phonon mean free paths is 
cut due to the presence of the heat reservoirs. This effect 
explains why finite size effects are less severe in NEMD 
than in EMD. We have also shown that the interfacial 
conductance measured in an EMD simulation whether 
using the Puech formula or the energy correlation func- 
tion identifies with the Landauer conductance which as- 
sumes phonons on both sides of the interface to be at 
equilibrium. This explains why in EMD a finite conduc- 
tance is measured when the two solids are similar. On 
the other hand, we have explained that in NEMD simula- 
tions, we measure a conductance given by the general ex- 
pression eq. (I10p inspired by Simons, and which accounts 
for the out-of-equilibrium distribution of phonons con- 
sistent with the imposed heat flux. Hence, we conclude 
that the two methods give intrinsically different values of 
the interfacial conductance. For impedance ratios typi- 
cal of real interfaces, the difference in conductances is 
large corresponding typically to a factor between 4 and 
10. On the other hand, when the impedance ratio is 
large the difference in conductances is small. This ex- 
plains why Barrat found good agreement between EMD 
and NEMD in the case of solid/liquid interfaces. Also we 
have shown that the two methods probe different energy 
transmission coefficients: EMD conductance are consis- 
tent with transmission describing diffuse events whose 
rates are governed primarily by the density of states mis- 
match between the two solids. On the other hand, in 
NEMD the transmission of phonons probed is specular 
in nature at least in the case analyzed here of atomically 
perfect interfaces. This difference stems from the differ- 
ent origin of the flux instantaneously flowing across the 
interface. 

An important question that we have to answer is which 
technique should be used-NEMD or EMD-to access 
a conductance measured experimentally. Intuitively, 
NEMD should be used to compute the value of the 
conductance measured experimentally using steady state 
technique such as the 3 omega method. On the other 
hand, EMD should be more akin to Laser pump-probe 
experiments where the transient response to an initial 
heating is recorded^. This needs further theoretical anal- 
ysis and will be the subject of future investigation. 
Another interesting question deals with the role of bal- 
listic phonons in the derivation of the non-equilibrium 



conductance eq. (fTOj) . Indeed, Landry and McGaughey 
observed that the non-equilibrium conductance overesti- 
mates the conductance measured at the interface between 
Si and Ge. We think that this discrepancy stems from the 
large value of the dominant mean free path in Si which 
is comparable with the system size considered. 
All these results have been obtained for the case of the 
perfect interfaces. This may allow to disentangle effects 
related to the contrast between the vibrational proper- 
ties of the bulk media from the effect arising from the 
interface. In particular, at the interface between real 
materials eventhough the interface may be treated to be- 
come atomically sharp, there is always a lattice mismatch 
which may enhance diffuse phonon scattering. The use of 
MD models allows then to measure each effect separately 
thus opening the way to a fundamental understanding of 
interfacial heat transfer between solids. 



VIII. APPENDIX: FINITE SIZE EFFECTS IN 
THE DETERMINATION OF THE EMD 
CONDUCTANCE 

In this appendix, we derive the length-dependent con- 
ductance eq.EHlmeasured in the EMD simulation. Please 
refer to the figure [§] for the relevant notations to be 
used here. As explained in the main body of the text, 
the length dependance of the conductance measured in 
EMD simulations is assumed to be caused by cross- 
correlations between the fluxes across the two interfaces 
of the system. We focus then on the cross correlation 
term (qx(t)qx' (0)) . We write the interfacial fluxes in 
terms of the phonons distribution functions in medium 
i: /^(r,*) = / e |(f) + Sf.^(f,t) where we have omitted 

the index designating the polarization to simplify the di- 
cusssion. Hence, we have : 

q A ,(t = 0)= i t2iSf 2 j;(0,r // ,t = Q)twv2xdf// 
+ / ^ Yl ti 2 5f lU (0,^ // ,t = 0)hujviW// (55) 

JA ' V «lx<0 

A similar equation holds for <?A(t) but in the following, 
we will use the following equation which derives from the 
continuity of the interfacial flux: 

A % 

= / ^X^A.fctV 2 ,^//^)^!^// (56) 
A fe 

The cross correlation (q^(t)qA' (0)) will thus involve cor- 
relation of the phonon distribution function of the form 
(5f it %(0,r,t = 0)Sf a ,(L/2,r ,£)). We assume the ther- 
mally induced phonon modes propagating in the medium 
i to be incoherent and characterized by a mean free path 
A. £. Under these conditions, the phonon correlation 
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writes: 



(W-t%(?>o)) 



Vit 



(57) 



where 



c v k B T 2 



2V(Huj) 2 

The cross correlation flux writes then: 



(« A '(*)?a(0)) 



EAksT 2 c v 2 { L 
— r t 2 iv 2x d [ — 



5 >o 



2V 



and the contribution to the conductance is 

c v t 2 ±v 2 

„>o 



(58) 

rj, /- ) exp(— |v2|i/A2) 
(59) 
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(60) 

To evaluate this latter conductance, we transform the 
discrete sum in a integral over the frequency: 

(q A ,{t)q A (0))dt 



FIG. 11. (Color online) Comparison between the exact Lan- 
dauer AMM conductance eq. [43] with the approximate solu- 
ex p (-L /2 cos 6A 2 )ti° n eq. (|44[) for the series of mass-mismatched Lennard- Jones 
crystals considered in the simulations. The mean acoustic ve- 
locity is taken to be that of the Lennard- Jones Argon. 
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Ak B T 2 Jo 
g 2 (u))c v t 21 {uj)\v 2 \I 2 (L,u)duj 



(61) 



where we have supposed that the transmission coefficient 
t 2 i is independent on the incidence angle and we have 
introduced the integral: 



I 2 (L,lu) 



exp 



Lu 



2A 2 (w) 



du 



For thick media, I » A and one can approximate the 
integral I(L,w) ~ exp (— ^) 45 - To evaluate the con- 
ductance eq. (l6l"j) . we assume as stated in the main body 
of the text that the vDOS is described by the Debye 
model and the frequency dependence of the mean free 
path is given by Callaway law22: 



A 2 (w) = A 2 \v 2 \/uj 2 



(63) 



where v 2 is assumed to be constant consistently with our 
hypothesis of Debye solid. The constant A 2 is related to 
the thermal conductivity X 2 through: 



A 2 = 



k B A 2 bj D ,2 
2n 2 \v 2 \ 



(64) 



If we suppose furthermore that the transmission coeffi- 
cient £12 is independent on the frequency uj as in the 
DMM model, it comes: 



Ak B T 2 



(QA'(t)qA(0))dt = -k B n 2 t 21 v 2 



2 1 ' Z A 2 



2/3 



TYV 2 L 



3/2 



where we have assumed that 0J max 3> \/v 2 A 2 j L, which 
physically means that the phonon mean free path of the 



mode with a frequency w max is smaller than the system 
length. This latter contribution may be rewritten: 



5G{L) 



3n 2 k B t 21 v 2 



A, 



G 22 L 



3/2 



(65) 



where we have introduced the conductance G 22 = 
3n 2 k B v 2 /8. Similar calculations allow to express the to- 
tal contribution of the cross fluxes as: 



(62) SG(L) 



n 2 t 2 iv 2 



\G 22 L 



3/2 



+ nit± 2 vi 



Ai 



3/2 



(66) 



Note that for Lennard- Jones solids differing only by their 
mass, the length £ = Xi/Gu = Ai/{n^ 3 Vi) is constant 
independent on the mass. Again anticipating the results 
of the section fVI El we can assume that the infinite length 
conductance Gi 2 is given by Gi 2 = \k B n 2 v 2 t 2 \ and the 
transmission coefficient obeys: t\ 2 = m 2 t 2 i/mi where we 
have introduced the masses of the two solids. Hence for 
the interface considered in fig. [5] for which m 2 /mi = 2, 
the correction to the conductance writes: 



5G{L) = cGi 2 f 



where c = J\ (1 + 2y/2) ~ 2.75. 



IX. 



3/2 



(67) 



APPENDIX: APPROXIMATIONS OF THE 
AMM CONDUCTANCES 



In this appendix, we assess the accuracy of different ap- 
proximations used to estimate the conductance appear- 
ing in eqs. (|4"3"|) and ([ST]) . More specifically, one needs 
to approximate the three geometrical integrals which de- 
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FIG. 12. (Color online) Comparison between the exact non- 
equilibrium conductance eq. (|5ip with the approximate solu- 
tion eq. 1)52) 1 for the series of mass-mismatched Lennard- Jones 
crystals considered in the simulations. The mean acoustic ve- 
locity is taken to be that of the Lennard- Jones Argon. 



pend on the Rayleigh transmission coefficient eq[ 



h = 
h = 

h = 



Mi*i2(j"i)dMi 



AtlM2*12(Ml)d/il 



(68) 
(69) 
(70) 



The general idea is to assume that the geometric inte- 
grals are dominated by phonons propagating in the soft 
material with a small angle 62 — when the acoustic 
mismatch betwen the two solids is large. This leads to 



use the following approximations: 
M2 — 1 

A*2 



1-21 

2 



M2 



1 



a 



(71) 
(72) 

(73) 



where we have denoted by a the ratio of the sound veloci- 
ties : a — C2/C1 < 1. The first approximation eq. (|7ip has 
been already discussed in the text and the correspond- 
ing approximate integrals are given in eqs. P5|) and (|53|) . 
in the second approximation eq. (|72[) . the three approxi- 
mated integrals depend on the parameter 7 = (3(1 — ^-): 



1 + 7 



+ 7 — 27 log 
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i-appx 
1 2 
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7 2 + 



(1+7) 



+ 37 log 



1 + 7 



7 appx = J ]_ . 



Q I 7-appx 

T^ 1 



4) 
(75) 



The third approximation eq. ()73j) yields calculations a 
little bit more involved. Within this approximation, one 
obtains the following expressions for the three integrals: 
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The denominator appearing in the three integrals has two 
poles 7*1/2 having multiplicity two and which are given 
by: 



(76) 
dp, (77) 
■dfi; (78) 
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(79) 
(80) 



and the approximated integral / appx are given by : 



/a ppx = 4/? / ^ 
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an 



hi 



4r, 
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n(i - n) f l- r i+ i 
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where i £ 1,2 and we have noted i + 1 = 1 + (i mod (2)), 
i.c, 1 + 1 = 2; 2 + 1 = 1. We can now study the accuracy 
of the previous approximations by comparing the exact 
expressions for the conductance eqs. (j43|) and (f5Tj) with 
the approximate equations involving the three approxi- 
mations discussed above. The comparison is shown in the 
figures (fTT|) and (fl2|) for the series of mass mismatched 
Lennard-Jones solids analyzed in the simulations. Strik- 
ingly the different approximations seem to work quite 
well over a broad range of acoustic impendance ratios. 
The first approximation eq fl7"Tj) slighlty underestimates 
the Landauer AMM conductance when the impedance 
ratio tends towards 1, but the two other approximations 
describe accurately the Landauer conductance for the 
whole range of ratio. As for the non-equilibrium con- 
ductance, the three approximations work quite well when 
the impedance ratio is smaller than 0.8. Above 0.8, the 
approximations eqs (|71jl and (|72| respectively overesti- 
mate and underestimate the conductance. In particular, 



eq. (|7Tj) predicts the divergence of the conductance at a 
value of the impedance ratio < 1. On the other hand 
the approximation eq (|73[) predicts accurately the final 
divergence of the conductance up to ratios ~ 0.9. None 
of the approximation presented predicts the divergence 
of the conductance when the impedance ratio tends to- 
wards 1, but in practice it is not common to work with 
such large ratios. 



ACKNOWLEDGMENTS 

Simulations have been run at the " Pole Scientifique de 
Modelisation Numerique" de Lyon using the LAMMPS 
open source package^. We acknowledge interesting dis- 
cussions with P. Chantrenne, T. Albaret, J.-Y. Duquesne 
and S. Volz. 



* samy.merabia@univ-lyonl.fr 

1 G. Chen, Phys. Rev. B 57 (1998) 14958 

2 E.T. Swartz and R.O. Pohl, Rev. Mod. Phys. 61 (1989) 
605 

3 D.G. Cahill, W.K. Ford, K.E. Goodson, G.D. Mahan, A. 
Majumdar, H.J. Maris, R. Merlin and S.R. Philpot, J. 
,4pp. Phys. 93(2003) 793 

4 R.J. Stoner and H.J. Maris, Phys. Rev. B 48 (1993) 16373 

5 Lyeo H.K. and D.G. Cahill, Phys. Rev. B 73 (2006) 144301 

6 V. Juve, M. Scardamaglia, P. Maioli, A. Crut, S. Merabia, 
L. Joly, N. Del Fatti and F. Vallee, Phys. Rev. B 80 (2009) 
195406 

7 P.E. Hopkins, P.M. Norris and R.J. Stevens, ASME J. Heat 
Transfer 130 (2008) 022401 

8 P.E. Hopkins and P.M. Norris, ASME J. Heat Transfer 
131 (2009) 022402 

9 DA. Young and H.J. Maris, Phys. Rev. B 40 (1989) 3685 
W. Zhang, T.S. Fisher and N. Mingo, AMSE J. Heat 

Transfer 129 (2007) 483 



11 E.S. Landry and A.J.H. McGaughey, Phys. Rev. B 80 
(2009) 165304 

12 R.J. Stevens, L.V. Zhigilei and P.M. Morris, Int. J. Heat 
Mass Transfer 50 (2007) 3977 

13 P.K. Schelling, S.R. Phillpot and P. Keblinski, Phys. Rev. 
B 65 (2002) 144306 

14 A. Rajabpour and S. Volz, J. Appl. Phys. 108 (2010) 
094324 

15 J.-L. Barrat and F. Chiaruttini, Mol. Phys. 101 (2003) 
1605 

16 A.J.H. McGaughey and J. Li, Proc. IMECE 2006 ASME 
conference (2006) 

17 S. Simons, J. Phys. C7 (1974) 4048 

18 E.T. Swartz and R.O. Pohl, Appl. Phys. Lett. 51 (1987) 
2200 

19 S. Pettersson and G.D. Mahan, Phys. Rev. B 42 (1990) 
7386 

20 G. Chen, App. Phys. Lett. 82 (2003) 991 

21 A comment is of interest here. Indeed, the interfacial quan- 
tity between two media should be evaluated at the nomi- 



18 



nal position of the Gibbs dividing surface, typically at half 
distance between two atomic layers. Hence, the tempera- 
ture profile should be extrapolated in each medium to the 
position of this surface. In this case, one measures a van- 
ishing temperature drop and an infinite conductance at the 
interface between two similar media. On the other hand, 
if one estimates a temperature jump between two adja- 
cent monolayers, one obtains an interfacial conductance 
on the order of A/a where A is the bulk conductivity of the 
material considered and a is the distance between nearest 
neighbouring atomic layers. This value of the conductance 
is typically two orders of magnitude larger than the max- 
imal conductance predicted by the Landauer equation in 
the classical limit nksv where n is the number density and 
v is the solid acoustic velocity. 

22 J. Callaway, Phys. Rev. 113 (1959) 1046 

23 W.A. Little, Can. J. Phys. 37 (1959) 334 

24 R. Landauer, Phil. Mag. 21 (1970) 863 

25 J.M. Ziman, Electrons and Phonons (Oxford University 
Press, New- York, 2001) 

26 Indeed this assumption is not necessary. Eq. [7] can be 
shown to be valid when the characteristic length of the 
problem (the distance between the thermostats) is larger 
than the phonon mean free path. 

27 J. A. Katerberg, C.L. Reynolds and A.C. Anderson, Phys. 
Rev. B 16 (1977) 673 

28 S. Aubry, C.J. Kimmer, A. Skye and P.K. Schelling, Phys. 
Rev. B 78 (2008) 064112 

29 Y. Chalopin, K. Esfarjani, A. Henry, S. Volz and G. Chen, 
Phys. Rev. B 85 (2012) 195302 

30 J. Stephenson, Physica A 117 (1983) 593 

31 The derivative of f(t) = exp(— |t|/rg) may be calculated 
using f(t) = H(t)exp(~t/Tj:) + H{-t) exp(t/rg) where 
H(t) is the Heaviside function. Using H'(t) = 5(t) one ar- 
rives at f'(t) = (sgn(t)/r s ) eatp(-|t|/rj) + 2tf(t). Note that 
indeed the derivative should be understand here in terms 



of generalized functions (distributions) 

32 This scaling is obtained because there are oc V terms in the 
sum over the allowded modes and because oc V/A. The 
thermodynamic limit in our case means A — > oo, V — > oo 
and V/A — > oo. As usual with Green-Kubo formulae, it 
is important to consider first the limit V — > oo before the 
limit t — > oo. 

33 K. Termentzidis, P. Chantrenne, J-Y. Duquesne and A. 
Saci, J. Phys. Cond. Matt. 22 (2010) 475001 

34 K. Termentzidis, J. Parasuraman, C.A. Da Cruz, S. Mer- 
abia, D. Angelescu, F. Marty, T. Bourouina, X. Kleber, P. 
Chantrenne and P. Basset, Nanoscale Research Letters, 6 
(2011) 288 

35 M.T. Dove, Introduction to Lattice Dynamics (Cambridge 
University Press, Cambridge England, 1993) 

36 N.W.Ashcroft and N.D. Mermin, Physique des Solides 
(EDP Sciences, Les Ulis, 2002) 

37 To arrive at eq. 1511 we have also used 

Jo Vl - (C2/ci)2 ^2i(^)dM2 = (ff) a /o A*i/*fl*«0*i)d/ii. 

38 L. Puech, G. Bonfait and B. Castaing, J. Low Temp. Phys. 
62 (1986) 315 

39 P. Chantrenne and J.-L. Barrat, J. Heat Transfer- 
Transactions of the ASME 126 (2004) 577 

40 D. Frenkel and B. Smit, Understanding Molecular simula- 
tion: from algorithms to applications Academic Press 2002 

41 P. Espanol and I. Zuniga, J. Chem. Phys. 98 (1993) 574 

42 K. Termentzidis, S. Merabia, P. Chantrenne and P. Ke- 
blinski, Int. J. Heat Mass Transf. 54 (2011) 2014 

43 A.J.H. McGaughey and M. Kaviany, Phys. Rev. B 69 
(2004) 094303 

44 K. Termentzidis, P. Chantrenne and P. Keblinski, Phys. 
Rev. B 79 (2009) 214307 

45 M.F. Modest, Radiative Heat Transfer^ Academic Press, 
Burlington USA, 2003) 

46 S. Plimpton, J. Comp. Phys. 117 (1995), 1-19: see 
http://lammps.sandia.gov 



